Most analyses conducted in this document are based on Capblancq and Forester (2021) and the associated Github repository.
Genomic data
# Population-based allele frequencies
# ===================================
geno <- read.csv(here("data/DryadRepo/ImputedGenomicData_AlleleFrequencies_withoutmaf.csv"),
row.names = 1)
# SNP sets
# ========
snp_sets <- readRDS(here("outputs/list_snp_sets.rds"))Climatic data
# Set of climatic variables
# =========================
clim_var <- readRDS(here("data/ClimaticData/NamesSelectedVariables.rds"))
# Climatic data
# =============
# we scale the past and future climatic data at the location of the populations
# with the parameters (mean and variance) of the past climatic data.
source(here("scripts/functions/generate_scaled_clim_datasets.R"))
clim_dfs <- generate_scaled_clim_datasets(clim_var,clim_ref_adj = FALSE)Projecting adaptive gradient(s) across space
Following Capblancq and Forester (2021) (and the associated Github repository), for the different sets of SNPs (one set of control SNPs and two sets of putative adaptive loci), a RDA was performed with the set of SNPs as multivariate response and the set of climatic variables as explanatory variables. For the sets of candidate SNPs, the resulting RDA spaces can be considered as ‘adaptively enriched genetic spaces’ (Capblancq and Forester 2021; Capblancq et al. 2023).
runRDA <- function(snp_set, clim_var, clim_df, geno){
# Keep the genomic data of the selected SNPs
geno <- geno[,snp_set]
# RDA formula
form_rda <- paste("geno ~ ", paste(clim_var,collapse= " + "))
# run RDA with only the selected SNPs
rda_outliers <- rda(as.formula(form_rda), clim_df)
}snp_sets <- lapply(snp_sets, function(x) {
x$rda_model <- runRDA(snp_set = x$set_snps,
clim_var = clim_var,
clim_df = clim_dfs$clim_ref,
geno = geno)
return(x)
}) %>% setNames(names(snp_sets))An RDA biplot can be used to visualize the relationship between the selected SNPs and the underlying climatic variables.
make_RDAbiplot <- function(x) {
TAB_loci <- as.data.frame(scores(x$rda_model, choices=c(1:2), display="species", scaling="none"))
TAB_var <- as.data.frame(scores(x$rda_model, choices=c(1:2), display="bp"))
eigenvalues <- summary(eigenvals(x$rda_model, model = "constrained")) %>% as.data.frame()
varexp_axis1 <- eigenvalues[2,1] *100
varexp_axis2 <- eigenvalues[2,2] *100
ggplot() +
geom_hline(yintercept=0, linetype="dashed", color = gray(.80), linewidth=0.6) +
geom_vline(xintercept=0, linetype="dashed", color = gray(.80), linewidth=0.6) +
geom_point(data = TAB_loci, aes(x=RDA1*3, y=RDA2*3), colour = '#CD24B2', size = 2, alpha = 0.8) + #"#F9A242FF"
geom_segment(data = TAB_var, aes(xend=RDA1, yend=RDA2, x=0, y=0), colour="black", linewidth=0.15, linetype=1, arrow=arrow(length = unit(0.02, "npc"))) +
geom_text(data = TAB_var, aes(x=1.1*RDA1, y=1.1*RDA2, label = row.names(TAB_var)), size = 3) +
xlab(paste0("RDA 1 (",round(varexp_axis1,1),"%)")) +
ylab(paste0("RDA 2 (",round(varexp_axis2,1),"%)")) +
facet_wrap(~paste0(x$set_name, " (n=",length(x$set_snps),")")) +
guides(color=guide_legend(title="Locus type")) +
theme_bw(base_size = 11) +
theme(panel.grid = element_blank(),
plot.background = element_blank(),
panel.background = element_blank(),
strip.text = element_text(size=11))
}rda_biplots <- lapply(snp_sets, make_RDAbiplot)
grid_RDAbiplots <- plot_grid(rda_biplots[[1]],rda_biplots[[2]],rda_biplots[[3]],
nrow=1)
ggsave(here("figs/RDA/RDAbiplots_SNPsets.pdf"), grid_RDAbiplots, width = 14, height=6)
# We save a plot for the Supplementary Information
# ================================================
plot_grid(rda_biplots[[1]],rda_biplots[[3]],nrow=1) %>% # we keep only the control and candidate SNPs
ggsave(here("figs/RDA/RDAbiplots_SNPsets_SI.pdf"), ., width = 9, height=6)
grid_RDAbiplots
We used the scores of the climatic variables along the RDA axes to calculate a genetic-based index for each pixel of the landscape. For each RDA axis, the index is calculated as follows:
\[\sum_{i=1}^{n}a_ib_i\]
where \(i\) is one of the \(n\) climatic variables used in the RDA model, \(a\) is the score of the climatic variable \(i\) along the RDA axis and \(b\) is the standardized value for the climatic variable \(i\) at the focal pixel.
When calculating based on the set of candidate SNPs, this index can be considered as an adaptive index that provides an estimate of the adaptive genetic similarity or difference of all pixels on the landscape as a function of the values of the climatic predictors at that location.
source(here("scripts/functions/project_adaptive_index.R"))
ai_proj <- lapply(snp_sets, function(snp_set){
project_adaptive_index(clim_var=clim_var,K=2,snp_set=snp_set)
})We project the index on maps to visualize the different gradients across the maritime pine range (which can be considered as adaptive gradients when based on the sets of putatively adaptive loci).
# countries borders for the map
world <- ne_countries(scale = "medium", returnclass = "sf")
ai_maps <- lapply(snp_sets, function(snp_set){
RDA_proj <- ai_proj[[snp_set$set_code]]
RDA_proj <- lapply(RDA_proj, function(x) rasterToPoints(x))
# The adaptive index is scaled between 0 and 1
for(i in 1:length(RDA_proj)){
RDA_proj[[i]][,3] <- (RDA_proj[[i]][,3]-min(RDA_proj[[i]][,3]))/(max(RDA_proj[[i]][,3])-min(RDA_proj[[i]][,3]))
}
# formatting the dataframes for ggplot
TAB_RDA <- as.data.frame(do.call(rbind, RDA_proj[1:2]))
colnames(TAB_RDA)[3] <- "value"
TAB_RDA$variable <- factor(c(rep("RDA1", nrow(RDA_proj[[1]])), rep("RDA2", nrow(RDA_proj[[2]]))), levels = c("RDA1","RDA2"))
# map options
point_size=2
x_limits = c(-10, 15)
y_limits = c(31, 52)
ggtitle <- snp_set$set_name
ggplot(data = TAB_RDA) +
geom_sf(data = world, fill="gray98") +
scale_x_continuous(limits = x_limits) +
scale_y_continuous(limits = y_limits) +
geom_raster(aes(x = x, y = y, fill = cut(value, breaks=seq(0, 1, length.out=11), include.lowest = T))) +
scale_fill_viridis_d(alpha = 0.8, direction = -1, option = "A") +
xlab("Longitude") + ylab("Latitude") +
guides(fill=guide_legend(title="Genetic index")) +
facet_grid(~ variable) +
ggtitle(ggtitle) +
theme_bw(base_size = 11) +
theme(panel.grid = element_blank(), plot.background = element_blank(), panel.background = element_blank(), strip.text = element_text(size=11))
})
# We save the maps for the Github repository
# ==========================================
pdf(here("figs/RDA/AdaptiveIndex_maps.pdf"), width=10,height=6)
lapply(ai_maps, function(x) x)
dev.off()
# We save the maps for the Supplementary Information
# ==================================================
# Candidate SNPs
pdf(here("figs/RDA/AdaptiveIndexMap_CandidateSNPs_SI.pdf"), width=9,height=4.9)
ai_maps[[1]] + theme(plot.title = element_blank())
dev.off()
# Control SNPs
pdf(here("figs/RDA/AdaptiveIndexMap_ControlSNPs_SI.pdf"), width=9,height=4.9)
ai_maps[[3]] + theme(plot.title = element_blank())
dev.off()
# Show AI maps
# ============
lapply(ai_maps, function(x) x)
Predicting the genomic offset of the studied populations under future climates.
# Function to calculate the RDA genetic offset for spatial points
# ===============================================================
# clim_ref = climatic data used to fit the RDA model (the climate-of-origin of the populations)
# clim_pred = climatic data used for predictions (so either the future climate of the populations,
# or climate of the common gardens or the NFI plots)
# weights = if NULL, the adaptive index is not weighted by the relative importance of the RDA axes in
# explaining the variance.
# K = number of RDA axes used to calculate the genomic offset
# snp_set = list with at least the RDA models
pred_GO_spatial_points <- function(snp_set,
K,
clim_ref,
clim_pred,
clim_var,
weights = NULL,
CG=F){
# weights based on the variance explained by the different axes
weights <- (snp_set$rda_model$CCA$eig/sum(snp_set$rda_model$CCA$eig))[1:K]
list_AI <- lapply(list(clim_ref,clim_pred), function(df){
lapply(1:K, function(i){
# Below we calculate the adaptive index for the axis i
# For that, we multiply the value of the variables at the location of the population
# by the loadings of the variables along the axis i
ai <- df %>%
dplyr::select(any_of(clim_var)) %>%
apply(1, function(x) sum(x * snp_set$rda_model$CCA$biplot[,i]))
if(!is.null(weights)){ai <- ai * weights[i]}
}) %>%
setNames(paste0("RDA", 1:length(.))) %>%
as.data.frame()
})
if(CG==F){
eucloffset <- unlist(lapply(1:nrow(list_AI[[1]]), function(x) dist(rbind(list_AI[[1]][x,], list_AI[[2]][x,]), method = "euclidean")))
} else if(CG==T){
eucloffset <- lapply(1:nrow(list_AI[[2]]), function(x) as.matrix(dist(rbind(list_AI[[2]][x,],list_AI[[1]]), method = "euclidean"))[2:(nrow(list_AI[[1]])+1),1]) %>%
setNames(clim_pred[,1] %>% pull()) %>%
as.data.frame() %>%
cbind(clim_ref[,1])
}
return(eucloffset)
}snp_sets <- lapply(snp_sets, function(x){
x$go <- lapply(names(clim_dfs$clim_pred), function(gcm){
pred_GO_spatial_points(snp_set = x,
clim_var=clim_var,
K = 2, # why K=2??
clim_ref = clim_dfs$clim_ref,
clim_pred = clim_dfs$clim_pred[[gcm]],
weights = T)
}) %>% setNames(names(clim_dfs$clim_pred))
return(x)
})source(here("scripts/functions/make_eucli_plot.R"))
# Calculate the Euclidean climatic distance
list_dist_env <- clim_dfs$clim_pred %>% lapply(function(clim_pred){
Delta = clim_dfs$clim_ref %>% dplyr::select(any_of(clim_var)) - clim_pred %>% dplyr::select(any_of(clim_var))
dist_env = sqrt( rowSums(Delta^2) )
})
# Main gene pools (for the figures)
gps <- readRDS(here("data/GenomicData/MainGenePoolPopulations.rds")) %>% arrange(pop)par(mfrow=c(1,2))
lapply(snp_sets, function(x) {
lapply(names(list_dist_env), function(gcm){
make_eucli_plot(
X = list_dist_env[[gcm]],
Y = x$go[[gcm]],
colors = gps$color_main_gp_pop,
color_names = gps$main_gp_pop,
ylab = "GDM genomic offset",
legend_position="topleft",
plot_title = paste0(x$set_name," - ", gcm))
})
}) 
# We generate scatter plots for the Supplementary Information.
# ============================================================
# Axis limits
# ===========
max_go <- lapply(snp_sets[c(1,3)], function(z){
z$go %>% unlist()
}) %>% unlist() %>% max()
max_go <- max_go + 0.01
range_eucli <- list_dist_env %>% unlist() %>% range()
# Run the function
# ================
lapply(snp_sets[c(1,3)], function(set_i) {
p <- lapply(names(list_dist_env), function(gcm){
make_ggscatterplot(
x = list_dist_env[[gcm]],
y = set_i$go[[gcm]],
title=gcm,
range_eucli = range_eucli,
max_go = max_go)
})
# remove y-labels to graphs in the second column
p[[2]] <- p[[2]] + ylab("")
p[[4]] <- p[[4]] + ylab("")
# remove x-labels to graphs in the second and third rows
p[[1]] <- p[[1]] + xlab("")
p[[2]] <- p[[2]] + xlab("")
p[[3]] <- p[[3]] + xlab("")
p[[6]] <- get_legend(p[[1]])
for(i in 1:5){p[[i]] <- p[[i]] + theme(legend.position = "none")}
plot_grid(plotlist=p, nrow = 3) %>%
ggsave(here(paste0("figs/RDA/ScatterPlotEucliDistance_",set_i$set_code,".pdf")),
.,
width=7,
height=8,
device="pdf")
})$cand
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/GOPredEvalPinpin/figs/RDA/ScatterPlotEucliDistance_cand.pdf"
$control
[1] "/home/juliette/Documents/Projects/GOPredEvalPinpin/GOPredEvalPinpin/figs/RDA/ScatterPlotEucliDistance_control.pdf"# Function to make the genomic offset maps
source(here("scripts/functions/make_go_map.R"))
# Population coordinates
pop_coord <- readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_noADJ.rds")))[[1]]$ref_means %>% dplyr::select(pop,longitude,latitude)
# Find minimum and maximum values of genomic offset for the maps
go_limits <- lapply(snp_sets, function(x) {
lapply(names(list_dist_env), function(gcm){
x$go[[gcm]]
}) %>% unlist()
}) %>% unlist() %>% range()
# The minimum GO value is very very small, almost zero, so we fix it to zero.
go_limits[[1]] <- 0
# Generate the maps for each set of SNPs and each GCM
lapply(snp_sets, function(x) {
go_maps <- lapply(names(list_dist_env), function(gcm){
make_go_map(dfcoord=pop_coord,
snp_set = x,
gcm=gcm,
ggtitle=gcm,
go_limits = go_limits,
point_size = 3)
})
legend_maps <- get_legend(go_maps[[1]])
go_maps <- lapply(go_maps, function(y) y + theme(legend.position = "none"))
go_maps$legend_maps <- legend_maps
go_maps <-plot_grid(plotlist=go_maps)
# =====================================================
# We save the figures for the Supplementary Information
# =====================================================
if(x$set_code=="cand"){
ggsave(here("figs/RDA/GOMaps_PopLocations_CandidateSNPs_SI.pdf"), go_maps, width=10,height=6, device="pdf")
} else if(x$set_code=="control"){
ggsave(here("figs/RDA/GOMaps_PopLocations_ControlSNPs_SI.pdf"), go_maps, width=10,height=6, device="pdf")
}
# =========
# Add title
# =========
title <- ggdraw() +
draw_label(
x$set_name,
fontface = 'bold',
x = 0,
hjust = 0
) +
theme(plot.margin = margin(0, 0, 0, 7))
# merge title and plots
plot_grid(
title, go_maps,
ncol = 1,
rel_heights = c(0.1, 1))
})
For each GCM, we attribute the value 1 to the top five populations with the highest genomic offset and we attribute the value 0 to the other populations. We then count the number of 1 for each population, which gives the table and map below:
source(here("scripts/functions/make_high_go_pop_maps.R"))
high_go_pops <- make_high_go_pop_maps(pop_coord=pop_coord,
list_go = snp_sets$cand$go,
ggtitle="RDA",
nb_id_pop = 5) # number of selected populations
saveRDS(high_go_pops, file = here("outputs/RDA/high_go_pops.rds"))
high_go_pops[[1]] %>% kable_mydf| pop | longitude | latitude | GFDL-ESM4 | IPSL-CM6A-LR | MPI-ESM1-2-HR | MRI-ESM2-0 | UKESM1-0-LL | sum_go |
|---|---|---|---|---|---|---|---|---|
| SEG | -8.45 | 42.82 | 0 | 0 | 0 | 0 | 0 | 0 |
| CAD | -6.42 | 43.54 | 0 | 0 | 0 | 0 | 0 | 0 |
| ARM | -6.46 | 43.30 | 0 | 0 | 0 | 0 | 0 | 0 |
| SAC | -8.36 | 42.12 | 0 | 0 | 0 | 0 | 0 | 0 |
| LAM | -6.22 | 43.56 | 0 | 0 | 0 | 0 | 0 | 0 |
| SIE | -6.49 | 43.53 | 0 | 0 | 0 | 0 | 0 | 0 |
| PUE | -6.63 | 43.55 | 0 | 0 | 0 | 0 | 0 | 0 |
| ALT | -6.49 | 43.28 | 0 | 0 | 0 | 0 | 0 | 0 |
| CAS | -6.98 | 43.50 | 0 | 0 | 0 | 0 | 0 | 0 |
| OLB | -0.62 | 40.17 | 0 | 0 | 0 | 0 | 0 | 0 |
| MIM | -1.30 | 44.13 | 0 | 0 | 0 | 0 | 0 | 0 |
| PET | -1.30 | 44.06 | 0 | 0 | 0 | 0 | 0 | 0 |
| PLE | -2.34 | 47.78 | 0 | 0 | 0 | 0 | 0 | 0 |
| HOU | -1.15 | 45.18 | 0 | 0 | 0 | 0 | 0 | 0 |
| STJ | -2.03 | 46.76 | 0 | 0 | 0 | 0 | 0 | 0 |
| SAL | -3.06 | 41.83 | 0 | 0 | 0 | 0 | 0 | 0 |
| VER | -1.09 | 45.55 | 0 | 0 | 0 | 0 | 0 | 0 |
| TAM | -5.02 | 33.60 | 0 | 0 | 0 | 0 | 0 | 0 |
| OLO | -1.83 | 46.57 | 0 | 0 | 0 | 0 | 0 | 0 |
| BAY | -2.88 | 41.52 | 0 | 0 | 0 | 0 | 0 | 0 |
| PIE | 9.04 | 41.97 | 0 | 0 | 0 | 0 | 0 | 0 |
| CAR | -4.28 | 41.17 | 0 | 0 | 0 | 0 | 0 | 0 |
| LEI | -8.96 | 39.78 | 0 | 0 | 0 | 1 | 0 | 1 |
| CEN | -4.49 | 40.28 | 0 | 0 | 1 | 0 | 0 | 1 |
| VAL | -4.31 | 40.52 | 0 | 0 | 1 | 0 | 0 | 1 |
| BON | -1.66 | 39.99 | 0 | 0 | 0 | 0 | 1 | 1 |
| ORI | -2.35 | 37.53 | 0 | 0 | 0 | 0 | 1 | 1 |
| QUA | -0.36 | 38.97 | 0 | 0 | 0 | 1 | 1 | 2 |
| ARN | -5.12 | 40.19 | 0 | 0 | 1 | 1 | 0 | 2 |
| COC | -4.50 | 41.25 | 1 | 1 | 0 | 0 | 0 | 2 |
| CUE | -4.48 | 41.37 | 1 | 1 | 0 | 0 | 0 | 2 |
| PIA | 9.46 | 42.02 | 1 | 1 | 0 | 0 | 0 | 2 |
| MAD | -5.23 | 35.18 | 1 | 1 | 1 | 1 | 1 | 5 |
| COM | -3.95 | 36.83 | 1 | 1 | 1 | 1 | 1 | 5 |
high_go_pops[[2]]
We look at the correlation across the different genomic offset predictions at the location of the populations, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.
par(mfrow=c(1,3))
lapply(names(snp_sets[[1]]$go), function(gcm){
lapply(snp_sets, function(x) x$go[[gcm]]) %>%
as_tibble() %>%
cor() %>%
corrplot(method = 'number',type = 'lower',
diag = FALSE,mar=c(0,0,2,0),
title=gcm,
number.cex=2,tl.cex=1.5)
})
# Function to calculate the RDA genetic offset using rasters
# ==========================================================
# Arguments
# =========
# snp_set = list with at least the RDA models
# clim_var = selected climatic variables
# K = number of RDA axes used to calculate the genomic offset
# gcm = Global Climate Model of the raster with future climatic data
# clim_ref_adj = TRUE or FALSE, specify whether the point estimate climatic data used to scale the rasters should be adjusted or not for elevation
# ref_period = the reference period used to calculate the adaptive index, can be 1901-1950 or 1960-1991
# method = `loadings` or `predict`
pred_GO_rasters <- function(snp_set,
clim_var,
K,
gcm,
range_buffer = shapefile(here('data/Mapping/PinpinDistriEUforgen_NFIplotsBuffer10km.shp')),
method = "loadings",
clim_ref_adj = FALSE,
ref_period = "ref_1901_1950"){
# Load point estimate climatic data of the reference period
if(clim_ref_adj==TRUE){adj <- "ADJ"} else {adj <- "noADJ"}
clim_ref_pe <- readRDS(here(paste0("data/ClimaticData/MaritimePinePops/ClimatePopulationLocationPointEstimates_ReferencePeriods_",adj,".rds")))[[ref_period]]
# Extract scaling parameters, i.e. mean and variance
scale_params <- lapply(clim_var, function(x){
vec_var <- clim_ref_pe$ref_means[,x] %>% pull()
list(mean = mean(vec_var), sd = sd(vec_var))
}) %>% setNames(clim_var)
# Rasters with climatic data for the reference period
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",clim_ref_pe$range[[1]],"-",clim_ref_pe$range[[2]],"_Extent-JulietteA/"))
tif_paths <- lapply(clim_var, function(x) paste0(path,"/",x,".tif"))
ref_rasts <- raster::stack(tif_paths)
# Raster with future climatic data
path <- here(paste0("data/ClimaticData/ClimateDTRasters/1km_",gcm,"_2041-2070_ssp370_Extent-JulietteA/"))
tif_paths <- lapply(clim_var, function(x) paste0(path,"/",x,".tif"))
fut_rasts <- raster::stack(tif_paths)
# checking that the CRS is the same for the buffer and the rasters
if(identical(crs(range_buffer),crs(ref_rasts))==FALSE){
stop(paste0("CRS of the range buffer is not the same as the CRS of the reference period raster."))}
if(identical(crs(range_buffer),crs(fut_rasts))==FALSE){
stop(paste0("CRS of the range buffer is not the same as the CRS of future climate rasters."))}
# Mask with the range if supplied
if(!is.null(range_buffer)){
ref_rasts <- mask(ref_rasts, range_buffer)
fut_rasts <- mask(fut_rasts, range_buffer)}
# extract coordinates and climatic values, and mean-center the climatic data
ref_vals <- as.data.frame(rasterToPoints(ref_rasts))
for(i in clim_var){ref_vals[,i] <- (ref_vals[,i] - scale_params[[i]]$mean) / scale_params[[i]]$sd}
fut_vals <- as.data.frame(rasterToPoints(fut_rasts))
for(i in clim_var){fut_vals[,i] <- (fut_vals[,i] - scale_params[[i]]$mean) / scale_params[[i]]$sd}
# Predict genomic offset for each pixel based on the loadings of the climatic variables
if(method == "loadings"){
ref_proj <- list()
fut_proj <- list()
go_proj <- list()
# Projection for each RDA axis
for(i in 1:K){
# Calculate adaptive index for the reference period
ref_rast <- ref_rasts[[1]]
ref_rast[!is.na(ref_rast)] <- as.vector(apply(ref_vals[,clim_var], 1, function(x) sum( x * snp_set$rda_model$CCA$biplot[,i])))
names(ref_rast) <- paste0("RDA_ref_", as.character(i))
ref_proj[[i]] <- ref_rast
names(ref_proj)[i] <- paste0("RDA", as.character(i))
# Calculate adaptive index under future climates
fut_rast <- fut_rasts[[1]]
fut_rast[!is.na(fut_rast)] <- as.vector(apply(fut_vals[,clim_var], 1, function(x) sum( x * snp_set$rda_model$CCA$biplot[,i])))
names(fut_rast) <- paste0("RDA_fut_", as.character(i))
fut_proj[[i]] <- fut_rast
names(fut_proj)[i] <- paste0("RDA", as.character(i))
# Calculate genetic offset based on a single RDA axis
go_proj[[i]] <- abs(ref_proj[[i]] - fut_proj[[i]])
names(go_proj)[i] <- paste0("RDA", as.character(i))
}
}
# Predict genomic offset for each pixel based on predict.RDA
if(method == "predict"){
# Prediction with the RDA model under reference-period climates and future climates
ref_pred <- predict(snp_set$rda_model, ref_vals[,-c(1,2)], type = "lc")
fut_pred <- predict(snp_set$rda_model, fut_vals[,-c(1,2)], type = "lc")
ref_proj <- list()
fut_proj <- list()
go_proj <- list()
for(i in 1:K){
# Calculate adaptive index for the reference period
ref_rast <- rasterFromXYZ(data.frame(ref_vals[,c(1,2)], Z = as.vector(ref_pred[,i])), crs = crs(ref_rasts))
names(ref_rast) <- paste0("RDA_ref_", as.character(i))
ref_proj[[i]] <- ref_rast
names(ref_proj)[i] <- paste0("RDA", as.character(i))
# Calculate adaptive index under future climates
fut_rast <- rasterFromXYZ(data.frame(ref_vals[,c(1,2)], Z = as.vector(fut_pred[,i])), crs = crs(ref_rasts))
names(fut_rast) <- paste0("RDA_fut_", as.character(i))
fut_proj[[i]] <- fut_rast
names(fut_proj)[i] <- paste0("RDA", as.character(i))
# Calculate genetic offset based on a single RDA axis
go_proj[[i]] <- abs(ref_proj[[i]] - fut_proj[[i]])
names(go_proj)[i] <- paste0("RDA", as.character(i))
}
}
# Weights based on axis eigen values
weights <- snp_set$rda_model$CCA$eig/sum(snp_set$rda_model$CCA$eig)
# Weight the current and future adaptive indices based on the eigen values of the associated axes
ref_proj_weighted <- do.call(cbind, lapply(1:K, function(x) rasterToPoints(ref_proj[[x]])[,-c(1,2)]))
ref_proj_weighted <- as.data.frame(do.call(cbind, lapply(1:K, function(x) ref_proj_weighted[,x]*weights[x])))
fut_proj_weighted <- do.call(cbind, lapply(1:K, function(x) rasterToPoints(fut_proj[[x]])[,-c(1,2)]))
fut_proj_weighted <- as.data.frame(do.call(cbind, lapply(1:K, function(x) fut_proj_weighted[,x]*weights[x])))
# Predict a global genetic offset, incorporating the K first axes weighted by their eigen values
go_global_proj <- go_proj[[1]]
go_global_proj[!is.na(go_global_proj)] <- unlist(lapply(1:nrow(ref_proj_weighted), function(x) dist(rbind(ref_proj_weighted[x,], fut_proj_weighted[x,]), method = "euclidean")))
names(go_global_proj) <- "global_offset"
# Return projections for current and future climates for each RDA axis, prediction of genetic offset for each RDA axis and a global genetic offset
return(list(ref_proj = ref_proj, fut_proj = fut_proj, go_proj = go_proj, go_global_proj = go_global_proj, weights = weights[1:K]))
}# For each snp set and each GCM, we generate a list of rasters with:
# the projection of AI for reference_period and future climates
# genomic offset predictions for each RDA axis
# genomic offset predictions incorporating the K first axes weighted by their eigen values
go_rasters <- snp_sets %>% lapply(function(snp_set){
lapply(names(clim_dfs$clim_pred), function(gcm){
pred_GO_rasters(snp_set,clim_var,K=2,gcm)}) %>%
setNames(names(clim_dfs$clim_pred))
})
go_rasters %>% saveRDS(file=here("outputs/RDA/go_rasters.rds"))go_rasters <- readRDS(file=here("outputs/RDA/go_rasters.rds"))
# Map options
# ===========
point_size=2
x_limits = c(-10, 15)
y_limits = c(31, 52)
# get maximum and minimum genomic offset values across the different GO predictions (for each snp set and each GCM)
# so that comparing the maps will be easier
go_limits <- names(snp_sets) %>% lapply(function(x){
lapply(names(clim_dfs$clim_pred), function(gcm){
rasterToPoints(go_rasters[[x]][[gcm]]$go_global_proj) %>% as_tibble() %>% pull(global_offset)
}) %>%
unlist()
}) %>%
unlist() %>%
range()
# The minimum GO value is very very small, almost zero, so we fix it to zero.
go_limits[[1]] <- 0
# Mapping
# =======
go_maps <- names(snp_sets) %>% lapply(function(x){
go_maps <- lapply(names(clim_dfs$clim_pred), function(gcm){
# extract genomic offset values
go_dfs <- rasterToPoints(go_rasters[[x]][[gcm]]$go_global_proj) %>% as_tibble()
ggplot(data = go_dfs) +
geom_sf(data = world, fill="gray98") +
scale_x_continuous(limits = x_limits) +
scale_y_continuous(limits = y_limits) +
geom_raster(aes(x = x, y = y, fill = global_offset), alpha = 1) +
scale_fill_gradient2(low="blue", mid= "yellow", high="red",
midpoint=(max(go_limits[[2]])-min(go_limits[[1]]))/2,
limits=go_limits,
name = "Genomic offset") +
xlab("Longitude") + ylab("Latitude") +
ggtitle(gcm) +
theme_bw(base_size = 11) +
theme(panel.grid = element_blank(),
plot.background = element_blank(),
panel.background = element_blank(),
#legend.position = c(0.85,0.15),
legend.box.background = element_rect(colour = "gray80"),
strip.text = element_text(size=11))
})
legend_maps <- get_legend(go_maps[[1]])
go_maps <- lapply(go_maps, function(y) y + theme(legend.position = "none"))
go_maps$legend_maps <- legend_maps
go_maps <-plot_grid(plotlist=go_maps)
# We save the maps for the Supplementary Information
# ==================================================
if(x=="cand"){
ggsave(here("figs/RDA/GOMap_CandidateSNPs_SI.pdf"), go_maps, width=11,height=7, device="pdf")
} else if(x=="control"){
ggsave(here("figs/RDA/GOMap_ControlSNPs_SI.pdf"), go_maps, width=11,height=7, device="pdf")
}
# Add title
# =========
title <- ggdraw() +
draw_label(
snp_sets[[x]]$set_name,
fontface = 'bold',
x = 0,
hjust = 0) +
theme(plot.margin = margin(0, 0, 0, 7))
# merge title and plots
plot_grid(title, go_maps, ncol = 1,rel_heights = c(0.1, 1))
})
# We save the maps for the Github repository
# ==========================================
pdf(here("figs/RDA/GOmaps.pdf"), width=14,height=9.5)
lapply(go_maps, function(x) x)
dev.off()
# Show GO maps
# ============
lapply(go_maps, function(x) x)
We check that the genomic offset predictions we obtained with the spatial points are highly correlated (ie similar) to those obtained with the rasters.
Comment When we extract the genomic offset values from the rasters at the location of the populations, we have missing values for some populations. Indeed, some populations are not included in the buffer of the species range that I used, based on the EUFORGEN distribution and 10km around the NFI plots. We may want to modify the buffer of the species range to include those populations!
# checking correlations
names(snp_sets) %>% lapply(function(x){
cor_go <- lapply(names(clim_dfs$clim_pred), function(gcm){
go_rast <- raster::extract(go_rasters[[x]][[gcm]]$go_global_proj, pop_coord[,c("longitude","latitude")])
cor(snp_sets[[x]]$go[[gcm]],go_rast,use= "complete.obs")
}) %>% unlist()
tibble(gcm=names(clim_dfs$clim_pred),
cor_go=cor_go)
}) %>%
setNames(names(snp_sets)) %>%
bind_rows(.id="snp_set") %>%
kable_mydf()| snp_set | gcm | cor_go |
|---|---|---|
| cand | GFDL-ESM4 | 0.99 |
| cand | IPSL-CM6A-LR | 1.00 |
| cand | MPI-ESM1-2-HR | 1.00 |
| cand | MRI-ESM2-0 | 1.00 |
| cand | UKESM1-0-LL | 1.00 |
| cand_corrected | GFDL-ESM4 | 0.99 |
| cand_corrected | IPSL-CM6A-LR | 1.00 |
| cand_corrected | MPI-ESM1-2-HR | 1.00 |
| cand_corrected | MRI-ESM2-0 | 1.00 |
| cand_corrected | UKESM1-0-LL | 0.99 |
| control | GFDL-ESM4 | 0.99 |
| control | IPSL-CM6A-LR | 1.00 |
| control | MPI-ESM1-2-HR | 1.00 |
| control | MRI-ESM2-0 | 1.00 |
| control | UKESM1-0-LL | 1.00 |
This is ok!
# Load the climatic data of the NFI plots.
nfi_clim <- readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))
# Keep only the climatic variables of interest and scale the climatic data
nfi_dfs <- generate_scaled_clim_datasets(clim_var, clim_ref = nfi_clim$clim_ref, clim_pred = nfi_clim$clim_survey)
# calculate the genomic offset for the NFI plots
snp_sets <- lapply(snp_sets, function(x){
x$go_nfi <- pred_GO_spatial_points(snp_set = x,
K = 2, # why K=2??
clim_var=clim_var,
clim_ref = nfi_dfs$clim_ref,
clim_pred = nfi_dfs$clim_pred,
weights = T)
return(x)
})
# checking missing data
# lapply(snp_sets, function(x) sum(is.na(x$go_nfi)))
# Find minimum and maximum values of genomic offset for the maps
go_limits <- lapply(snp_sets, function(snp_set) snp_set$go_nfi) %>% unlist() %>% range()
# The minimum GO value is very very small, almost zero, so we fix it to zero.
go_limits[[1]] <- 0
# map genomic offset predictions in the NFI plots
lapply(snp_sets, function(x) {
p <- make_go_map(
dfcoord= readRDS(here("data/ClimaticData/NFIplots/NFIclimate.rds"))[[1]] %>% dplyr::select(contains("ude")),
snp_set = x,
type="NFI",
point_size = 0.5,
go_limits = go_limits,
legend_position = c(0.85,0.15),
legend_box_background = "gray80",
y_limits = c(35, 51))
# Figure for the SI
# =================
if(x$set_code=="cand"){
p_SI <- p + theme(plot.title = element_blank())
ggsave(here("figs/RDA/NFI_GOmap_CandidateSNPs_SI.pdf"), p_SI, width=8,height=8, device="pdf")
} else if(x$set_code=="control"){
p_SI <- p + theme(plot.title = element_blank())
ggsave(here("figs/RDA/NFI_GOmap_ControlSNPs_SI.pdf"), p_SI, width=8,height=8, device="pdf")
}
# Show maps in the Quarto document
# ================================
p
})
We look at the correlation across the different genomic offset predictions in the NFI plots, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.
lapply(snp_sets, function(x) x$go_nfi) %>%
as_tibble() %>%
cor() %>%
corrplot(method = 'number',type = 'lower', diag = FALSE,mar=c(0,0,2,0),
number.cex=2,tl.cex=1.5)
cg_clim <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% dplyr::select(cg,any_of(clim_var))
cg_coord <- readRDS(here("data/ClimaticData/CommonGardens/ClimateCG.rds")) %>% dplyr::select(cg,contains("ude"))
cg_names <- unique(cg_coord$cg)
# Generate scaled climatic datasets with climatic data at the location of the populations and at the location of the common gardens
cg_dfs <- generate_scaled_clim_datasets(clim_var, clim_pred = cg_clim)
# Predict genomic offset of each population when transplanted in the climate of the common gardens
snp_sets <- lapply(snp_sets, function(x){
x$go_cg <- pred_GO_spatial_points(snp_set = x,
K = 2, # why K=2??
clim_var = clim_var,
clim_ref = cg_dfs$clim_ref,
clim_pred = cg_dfs$clim_pred,
weights = T,
CG=T) %>% as_tibble()
return(x)
})
# Map genomic offset predictions at the locations of the populations
go_maps_cg <- lapply(cg_names, function(cg_name){
# Find minimum and maximum values of genomic offset for the maps
go_limits <- lapply(snp_sets, function(snp_set) snp_set$go_cg[[cg_name]]) %>% unlist() %>% range()
go_limits[[1]] <- 0
p <- lapply(snp_sets, function(x) make_go_map(dfcoord=pop_coord,
snp_set = x,
point_size = 3,
type="CG",
go_limits = go_limits,
cg_name=cg_name,
cg_coord=cg_coord))
plot_grid(p[[1]],p[[2]],p[[3]],nrow=1)
}) %>% setNames(cg_names)
pdf(here("figs/RDA/GOmaps_CGs.pdf"), width=15,height=6)
lapply(go_maps_cg, function(x) x)
dev.off()
# show maps
lapply(go_maps_cg, function(x) x)
We look at the correlation across the different genomic offset predictions in the common gardens, i.e. those based on all SNPs and those based on sets of candidates or control SNPs.
par(mfrow=c(1,3))
lapply(cg_names, function(cg_name){
lapply(snp_sets, function(x) x$go_cg) %>%
lapply(function(set){
set[[cg_name]]
}) %>%
as_tibble() %>%
cor() %>%
corrplot(method = 'number',type = 'lower',
diag = FALSE,mar=c(0,0,2,0),
title=str_to_title(cg_name),
number.cex=2,tl.cex=1.5)
})
Let’s save the genomic offset predictions for comparison with the other methods.
snp_sets %>% saveRDS(file=here("outputs/RDA/go_predictions.rds"))devtools::session_info()─ Session info ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
setting value
version R version 4.3.2 (2023-10-31)
os Ubuntu 22.04.3 LTS
system x86_64, linux-gnu
ui X11
language (EN)
collate fr_FR.UTF-8
ctype fr_FR.UTF-8
tz Europe/Paris
date 2024-01-18
pandoc 3.1.1 @ /usr/lib/rstudio/resources/app/bin/quarto/bin/tools/ (via rmarkdown)
─ Packages ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
package * version date (UTC) lib source
cachem 1.0.8 2023-05-01 [1] CRAN (R 4.3.2)
class 7.3-22 2023-05-03 [4] CRAN (R 4.3.1)
classInt 0.4-10 2023-09-05 [1] CRAN (R 4.3.2)
cli 3.6.2 2023-12-11 [1] CRAN (R 4.3.2)
cluster 2.1.6 2023-12-01 [4] CRAN (R 4.3.2)
codetools 0.2-19 2023-02-01 [4] CRAN (R 4.2.2)
colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.3.2)
corrplot * 0.92 2021-11-18 [1] CRAN (R 4.3.2)
cowplot * 1.1.2 2023-12-15 [1] CRAN (R 4.3.2)
DBI 1.2.0 2023-12-21 [1] CRAN (R 4.3.2)
devtools 2.4.5 2022-10-11 [1] CRAN (R 4.3.2)
digest 0.6.33 2023-07-07 [1] CRAN (R 4.3.2)
dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.3.2)
e1071 1.7-14 2023-12-06 [1] CRAN (R 4.3.2)
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.3.2)
evaluate 0.23 2023-11-01 [1] CRAN (R 4.3.2)
fansi 1.0.6 2023-12-08 [1] CRAN (R 4.3.2)
farver 2.1.1 2022-07-06 [1] CRAN (R 4.3.2)
fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.3.2)
forcats * 1.0.0 2023-01-29 [1] CRAN (R 4.3.2)
fs 1.6.3 2023-07-20 [1] CRAN (R 4.3.2)
generics 0.1.3 2022-07-05 [1] CRAN (R 4.3.2)
ggplot2 * 3.4.4 2023-10-12 [1] CRAN (R 4.3.2)
glue 1.6.2 2022-02-24 [1] CRAN (R 4.3.2)
gtable 0.3.4 2023-08-21 [1] CRAN (R 4.3.2)
here * 1.0.1 2020-12-13 [1] CRAN (R 4.3.2)
highr 0.10 2022-12-22 [1] CRAN (R 4.3.2)
hms 1.1.3 2023-03-21 [1] CRAN (R 4.3.2)
htmltools 0.5.7 2023-11-03 [1] CRAN (R 4.3.2)
htmlwidgets 1.6.4 2023-12-06 [1] CRAN (R 4.3.2)
httpuv 1.6.13 2023-12-06 [1] CRAN (R 4.3.2)
httr 1.4.7 2023-08-15 [1] CRAN (R 4.3.2)
jsonlite 1.8.8 2023-12-04 [1] CRAN (R 4.3.2)
kableExtra * 1.3.4 2021-02-20 [1] CRAN (R 4.3.2)
KernSmooth 2.23-22 2023-07-10 [4] CRAN (R 4.3.1)
knitr * 1.45 2023-10-30 [1] CRAN (R 4.3.2)
labeling 0.4.3 2023-08-29 [1] CRAN (R 4.3.2)
later 1.3.2 2023-12-06 [1] CRAN (R 4.3.2)
latex2exp * 0.9.6 2022-11-28 [1] CRAN (R 4.3.2)
lattice * 0.22-5 2023-10-24 [4] CRAN (R 4.3.1)
lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.3.2)
lubridate * 1.9.3 2023-09-27 [1] CRAN (R 4.3.2)
magrittr * 2.0.3 2022-03-30 [1] CRAN (R 4.3.2)
MASS 7.3-60 2023-05-04 [4] CRAN (R 4.3.1)
Matrix 1.6-3 2023-11-14 [4] CRAN (R 4.3.2)
memoise 2.0.1 2021-11-26 [1] CRAN (R 4.3.2)
mgcv 1.9-1 2023-12-21 [4] CRAN (R 4.3.2)
mime 0.12 2021-09-28 [1] CRAN (R 4.3.2)
miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.3.2)
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.3.2)
nlme 3.1-163 2023-08-09 [4] CRAN (R 4.3.1)
permute * 0.9-7 2022-01-27 [1] CRAN (R 4.3.2)
pillar 1.9.0 2023-03-22 [1] CRAN (R 4.3.2)
pkgbuild 1.4.3 2023-12-10 [1] CRAN (R 4.3.2)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.3.2)
pkgload 1.3.3 2023-09-22 [1] CRAN (R 4.3.2)
profvis 0.3.8 2023-05-02 [1] CRAN (R 4.3.2)
promises 1.2.1 2023-08-10 [1] CRAN (R 4.3.2)
proxy 0.4-27 2022-06-09 [1] CRAN (R 4.3.2)
purrr * 1.0.2 2023-08-10 [1] CRAN (R 4.3.2)
R6 2.5.1 2021-08-19 [1] CRAN (R 4.3.2)
ragg 1.2.7 2023-12-11 [1] CRAN (R 4.3.2)
raster * 3.6-26 2023-10-14 [1] CRAN (R 4.3.2)
RColorBrewer * 1.1-3 2022-04-03 [1] CRAN (R 4.3.2)
Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.3.2)
readr * 2.1.4 2023-02-10 [1] CRAN (R 4.3.2)
remotes 2.4.2.1 2023-07-18 [1] CRAN (R 4.3.2)
rlang 1.1.2 2023-11-04 [1] CRAN (R 4.3.2)
rmarkdown 2.25 2023-09-18 [1] CRAN (R 4.3.2)
rnaturalearth * 1.0.1 2023-12-15 [1] CRAN (R 4.3.2)
rnaturalearthdata * 0.1.0 2017-02-21 [1] CRAN (R 4.3.2)
rprojroot 2.0.4 2023-11-05 [1] CRAN (R 4.3.2)
rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.3.2)
rvest 1.0.3 2022-08-19 [1] CRAN (R 4.3.2)
scales 1.3.0 2023-11-28 [1] CRAN (R 4.3.2)
sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.3.2)
sf * 1.0-15 2023-12-18 [1] CRAN (R 4.3.2)
shiny 1.8.0 2023-11-17 [1] CRAN (R 4.3.2)
sp * 2.1-2 2023-11-26 [1] CRAN (R 4.3.2)
stringi 1.8.3 2023-12-11 [1] CRAN (R 4.3.2)
stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.3.2)
svglite 2.1.3 2023-12-08 [1] CRAN (R 4.3.2)
systemfonts 1.0.5 2023-10-09 [1] CRAN (R 4.3.2)
terra 1.7-65 2023-12-15 [1] CRAN (R 4.3.2)
textshaping 0.3.7 2023-10-09 [1] CRAN (R 4.3.2)
tibble * 3.2.1 2023-03-20 [1] CRAN (R 4.3.2)
tidyr * 1.3.0 2023-01-24 [1] CRAN (R 4.3.2)
tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.3.2)
tidyverse * 2.0.0 2023-02-22 [1] CRAN (R 4.3.2)
timechange 0.2.0 2023-01-11 [1] CRAN (R 4.3.2)
tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.3.2)
units 0.8-5 2023-11-28 [1] CRAN (R 4.3.2)
urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.3.2)
usethis 2.2.2 2023-07-06 [1] CRAN (R 4.3.2)
utf8 1.2.4 2023-10-22 [1] CRAN (R 4.3.2)
vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.3.2)
vegan * 2.6-4 2022-10-11 [1] CRAN (R 4.3.2)
viridisLite 0.4.2 2023-05-02 [1] CRAN (R 4.3.2)
webshot 0.5.5 2023-06-26 [1] CRAN (R 4.3.2)
withr 2.5.2 2023-10-30 [1] CRAN (R 4.3.2)
xfun 0.41 2023-11-01 [1] CRAN (R 4.3.2)
xml2 1.3.6 2023-12-04 [1] CRAN (R 4.3.2)
xtable 1.8-4 2019-04-21 [1] CRAN (R 4.3.2)
yaml 2.3.8 2023-12-11 [1] CRAN (R 4.3.2)
[1] /home/juliette/R/x86_64-pc-linux-gnu-library/4.3
[2] /usr/local/lib/R/site-library
[3] /usr/lib/R/site-library
[4] /usr/lib/R/library
──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────